knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(spatstat)
library(maptools) 
# library(sp)
library(raster) ### BEFORE tidyverse! b/c select()
library(tidyverse)
library(here)
library(sf)
library(tmap)
library(gstat)
library(stars)
library(janitor)

1. Introduction

This report will carry out a spatial analysis using tmap to create an interactive map showing the location of soil spill events in california using the CA DFW Oil Spill Incident Tracking database.

This report will also produce the creation of a finalized static choropleth map using ggplot in which the fill color for each county depends on the count of inland oil spill events by county for the 2008 oil spill data.

Data Source: The data used for analysis is from the CA DFW Oil Spill Incident Tracking dataset. The database system is designed to provide OSPR with quantified statistical data on oil spill response by OSPR field responders. The OSPR Incident Tracking Database System project was initiated to provide OSPR with oil spill incident data for statistical evaluation and justification for program planning, drills and exercise training and development, legislative analysis, budget preparation, to inform and educate the public and analyze OSPR’s overall spill preparedness and response performance. An “incident”, for purposes of this database, is “a discharge or threatened discharge of petroleum or other deleterious material into the waters of the state.”

Data can also be downloaded from the State of California Geoportal

# Reading in the data
oil <- read_sf(here("data", "Oil_Spill_Incident_Tracking_[ds394]", "Oil_Spill_Incident_Tracking_[ds394].shp"))

# Check the projection
st_crs(oil) 
## Coordinate Reference System:
##   User input: WGS 84 / Pseudo-Mercator 
##   wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["Popular Visualisation Pseudo-Mercator",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Web mapping and visualisation."],
##         AREA["World between 85.06°S and 85.06°N."],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]
# Read in the CA county data (TIGER shapefile):
ca_counties_sf <- read_sf(here("data/ca_counties"), layer = "CA_Counties_TIGER2016") %>% 
  janitor::clean_names() %>% 
  select(name)

# Check the projection
st_crs(ca_counties_sf) 
## Coordinate Reference System:
##   User input: WGS 84 / Pseudo-Mercator 
##   wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["Popular Visualisation Pseudo-Mercator",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Web mapping and visualisation."],
##         AREA["World between 85.06°S and 85.06°N."],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]

2. Plotting Oil Data on CA Counties

#testing a plot to make sure the plot runs
#ggplot() +
 # geom_sf(data = ca_counties_sf) +
  #geom_sf(data = oil)
# TSetting Tmap to plot:
tmap_mode("view")

#Mapping with the oil data and the county data
tm_shape(ca_counties_sf) +
  tm_polygons(alpha = 0) +
  tm_shape(oil) +
  tm_dots(col = "navyblue")

Figure 1: Figure 1 illustrates the distribution of oil spill events indicated in dark blue from the CA DFW Oil Spill Incident Tracking dataset across California. California county boundaries are indicated in black

3. Cholopleth Map of Inland Oil Spills by CA County

Wangling data

oil_inland <- oil %>% 
  filter(INLANDMARI == "Inland") #filtering to keep inland

oil_inland_county <- oil_inland %>% 
  st_join(ca_counties_sf)

oil_counts <- oil_inland_county %>% 
  group_by(name) %>% 
  summarize(total_spills = sum(!is.na(DFGCONTROL)))

Oil Spill Events by County in 2008

ggplot() +
  geom_sf(data = oil_counts, aes(color = total_spills))

X. Data Citation

Data source: Mark Lampinen, Department of Fish and Game, Office of Spill Prevention and Response, 916 322-4777 CA DFW Oil Spill Incident Tracking dataset

END TASK